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NUMERICAL SOLUTION OF SUPERSONIC THREE-DIMENSIONAL 
FREE-MIXING FLOWS USING THE PARABOLIC-ELLIPTIC 
NAVIER-STOKES EQUATIONS 

Richard S. Hirsh 
Langley Research Center 

SUMMARY 

A numerical method is presented for solving the parabolic -elliptic Navier-Stokes 
equations. The solution procedure is applied to three-dimensional supersonic laminar 
jet flow issuing parallel with a supersonic free stream. A coordinate transformation is 
introduced which maps the boundaries at infinity into a finite computational domain in 
order to eliminate difficulties associated with the imposition of free-stream boundary 
conditions. 

Results are presented for an approximate circular jet, a square jet, varying aspect 
ratio rectangular jets, and interacting square jets. The solution behavior varies from 
axisymmetric to nearly two-dimensional in character. For cases where comparisons of 
the present results with those obtained from shear layer calculations could be made, 
agreement was good. 


INTRODUCTION 

Free jet and wake flows are common aerodynamic phenomena. These flows are 
generally turbulent, and the calculation of two-dimensional or axisymmetric turbulent free 
jets or wakes is difficult (ref. 1) because of problems associated with turbulence modeling; 
higher order modeling (Reynolds stress equations) is necessary in many cases. For 
three-dimensional flows with two essential cross-plane velocities, very few calculations 
have been made. To assess the modeling procedures for three-dimensional flows, it is 
necessary to develop a numerical procedure for laminar flows, preferably in primitive 
variables, to allow ease of incorporation of the turbulence models. The purpose of this 
paper is to present such a procedure and demonstrate its applicability. 

The calculation of free-mixing flows has, in the past, been accomplished through use 
of the boundary-layer assumptions in the two-dimensional or axisymmetric Navier-Stokes 
equations. The accuracy and validity of these procedures have been well documented in 
the literature. (See refs. 1 and 2.) However, there are numerous situations where the 
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flow cannot be considered either two-dimensional or axisymmetric. Jets issuing from 
rectangular orifices, wakes behind any but the simplest bodies, and the flow downstream 
of a wing tip are examples of three-dimensional free-mixing flows where boundary-layer 
assumptions are invalid. The characteristic feature of these flows is the importance of 
diffusion in two spatial coordinates. 

These flows have certain characteristics which are common to boundary -layer flows; 
for example, the velocities in the plane normal to the main-stream direction are usually 
much smaller than the main-flow velocity. Also, one expects gradients that exist in the 
cross directions to be larger than the gradients in the main-flow direction. Therefore, it 
might be reasonable to use a boundary-layer scaling on the Navier-Stokes equations to 
effect some simplification; however, this procedure yields an inconsistent set of equations. 
To lowest order, the cross-stream momentum equations reduce to statements that the 
pressure is constant, and the resulting number of equations is not sufficient to determine 
the remaining unknown quantities. A set of equations which is a higher order approxi- 
mation to the Navier-Stokes than the boundary -layer equations, the parabolic -elliptic 
Navier-Stokes (PENS) equations, is necessary to overcome this difficulty. 

The PENS equations were originally derived (in two dimensions) in reference 3 for 
the merged flow region near the leading edge of a flat plate in hypersonic flow. An explicit 
finite -difference method was used to integrate the equations in the domain between the flat 
plate and the free -stream flow beyond the merged shock layer. The three-dimensional 
equations are given in reference 4 where the leading-edge flow on a finite -width flat plate 
is calculated. The alternating direction implicit (ADI) method, introduced in reference 5, 
is used to solve the equations. The two- and three-dimensional problems are reviewed in 
reference 6 where a new predictor -corrector method is presented as the most appropriate 
one to use. See also reference 7 for further details. This analysis has been extended to 
the hypersonic tip region of a cone in reference 8. 

The flow past a cone received further study in references 9 and 10 where the PENS 
equations were used to calculate the entire flow field between a well defined shock envelope 
and the cone surface. The equations solved in references 9 and 10 are very similar to 
those of reference 8, but the arguments leading to their derivation are somewhat different. 
A fully implicit numerical method is used to solve the resulting equations between the body 
and the shock, where the Rankine-Hugoniot equations are imposed. 

The PENS concepts have also been utilized in subsonic, incompressible flows for 
internal (ref. 11) and external (ref. 12) cases. However, the need to solve a Poisson 
equation for the cross-plane pressure in these incompressible flows separates them from 
the procedures used in references 3, 4, and 6 to 10 and the present work where only super- 
sonic flows are considered. 
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The present work solves the parabolic -elliptic Navier-Stokes equations for super- 
sonic free-mixing flows. It is found necessary to incorporate a coordinate transformation 
into the PENS equations in order to impose free-stream boundary conditions on the flow. 
The character of the equations is such that they are well suited to solution by the ADI 
method. This finite -difference method generates block tridiagonal matrix equations which 
allow for the simultaneous solution of all the variables at each location. Preliminary 
results from this investigation were given in reference 13. Further test cases and dis- 
cussion of results are presented in this paper. 

This solution procedure of the PENS system is applied to three-dimensional, super- 
sonic, rectangular jet flows in which the jet aspect ratio is varied from one (square jet) to 
large values representative of slits (two-dimensional jets). Various combinations of free- 
stream Mach and Reynolds numbers are studied. Whenever possible, results from free 
shear layer calculations obtained by using the boundary-layer approximations are com- 
pared with those obtained by using the PENS equations. The range of applicability of the 
procedure is demonstrated from the near axisymmetry of the square, through a true 
three-dimensional region of moderate aspect ratios, to a quasi-two-dimensional flow in 
the high-aspect-ratio limit which approximates a two-dimensional jet. 

Finally, a test case corresponding to interacting jet flow behind an aircraft was 
computed. No difficulties were encountered during the course of these calculations which 
gave the correct qualitative results. No quantitative comparisons could be made with 
any two-dimensional or axisymmetric shear-layer predictions since this flow is three- 
dimensional, and beyond the scope of boundary -layer theory. 

SYMBOLS 


A,B transformation constants 

[a], [b], [c] matrices defined by equation (23) 

a,b,c components of matrices [a], [b], [c] 

a,b,a coefficients in equation (Al) 

Cp specific heat 

D vector defined by equation (23) 

— convective derivative 

Dt 
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d component of vector D 

h width of unit jet (reference length) 

L r reference length 

M Mach number 

Np r Prandtl number 

Np e Reynolds number 

p pressure 

Q any representative dependent variable 

5 =S*/T r 

S* Sutherland's constant 

T temperature 

u,v,w component velocities in x-, y-, and z -directions, respectively 

W vector of dependent variables (see eq. (24)) 

x,y,z coordinate directions 

y ratio of specific heats 

A dilatation 

Ax,Ay,Az grid spacings in x-, y-, and z-directions, respectively 
At),A£ grid spacings in?]- and ^-directions, respectively 

6 central difference operator; also measure of boundary-layer scale 

7],£ transformed y- and z -coordinates 
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M 


viscosity 


£ streamwise coordinate in transformed plane 

p density 

$ dissipation 

Superscripts: 

i discretized x-position 

* intermediate step of alternating direction implicit method 

T transpose 

Subscripts: 

j,k discretized y- and z -positions, respectively 

r reference quantity 

°° free stream 

+, - upper and lower grid spacings used with nonuniform grid 

<t. center line 


A bar over a symbol denotes a nondimensional quantity. 

GOVERNING EQUATIONS 


The Navier -Stokes equations for the steady, three-dimensional flow of a viscous, 
heat-conducting fluid are, in Cartesian coordinates 

for x- momentum, 
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for z- momentum, 
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and the continuity equation is 


4 (p u) + ^ {PV) + ^ (PW) =0 (5) 

The elliptic nature of these equations makes it necessary, for a numerical solution, 
to store the unknowns for the entire three-dimensional region being computed. The core 
storage available on present computers is insufficient to practicably handle any but the 
coarsest grids, which yield computed results of very low accuracy. Thus, methods to 
reduce the equations to a form more tractable for computations must be employed. 


Boundary-Layer Approximation 

The obvious simplification to be tried, for viscous flows with a dominant flow 
direction, is the classical boundary-layer technique which scales the various field quan- 
tities and their gradients according to their relative sizes. This approximation will be 
seen to be invalid for three-dimensional flows where, as assumed below, the cross-plane 
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velocities and characteristic length scales are of the same order of magnitude, but not 
necessarily equal. 

If the x-direction is assumed to be the dominant direction of flow, and the main 
changes in the flow field are taken to occur perpendicular to this main stream, that is, 
in the YZ plane, the boundary-layer scaling for the equations is 
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Placing this scaling into equations (1) to (4) gives (after dropping the overbar) 
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where the energy equation has been omitted since it adds nothing to the following argument 
concerning the momentum equations. Making the boundary-layer approximation by 
retaining only the highest order terms in the Reynolds number gives 
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4 (pu) + ^ (pv) + 4 (pw) = 0 


Hence, the two cross momentum equations have degenerated into a statement that the 
pressure is constant in cross planes and the system of equations which results is not 
sufficient to determine all the unknowns. This condition was noted previously by Greber 
in a comment following an article by Pai and Hsieh (ref. 14) who used an a priori boundary - 
layer scaling on equations similar to equations (1) to (3). Hence, a true boundary-layer 
scaling cannot be expected to yield one simple set of equations which can be easily solved 
and a different approximation must be made. 


Parabolic -Elliptic Navier- Stokes Equations 


Although a strict application of boundary-layer scaling is abandoned, some of the 
concepts from the boundary-layer theory indicate a means to simplify the equations. The 
only assumption which can be made is the predominance of the convection in the one main- 
flow direction over the shear stress in this direction. This assumption leads to the 
(Reynolds number dependent) conclusion that x-derivatives associated with the shear 
stress can be neglected. With the x-direction taken to be the convective main-stream 
direction, the Navier -Stokes equations, after neglecting all x-derivatives associated with 
shear, become 
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The x-momentum equation is the same as would have been produced by a true 
boundary -layer scaling, but since no quantitative assumptions have been made concerning 
the relative sizes of the x-, y-, or z -derivatives, the y- and z -momentum equations are 
retained, although in somewhat simpler form than equations (2) and (3). The continuity 
equation remains unchanged, that is, equation (5). With the same arguments, the energy 
equation becomes 
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When these equations are supplemented by an equation of state and a viscosity relation 


p = PRT 

4 far + S*\ /T \ 3 / 2 
yT + S* AtJ 


where R is the gas constant, a system of five equations for five unknowns is obtained 
after elimination of the density by the perfect gas equation of state. 

The elliptic nature of the Navier-Stokes equations in the x-direction has thus been 
eliminated; consequently, the equations are parabolic in x and marching integration may 
be used in the streamwise direction. This simplification significantly reduces the 
computational requirements since it eliminates the need to store all the field quantities 
at each x-location and results in a substantial reduction in computer storage. Thus, 
they are given the name parabolic -elliptic Navier-Stokes equations, since the assumptions 
allow a march in x away from an initial data plane and yet retain the elliptic character 
of the crossflow planes (YZ planes) due to the inclusion of all second derivatives in 
y and z. For example, flows with swirl or possible crossflow recirculation (vortices) 
in the YZ planes can be computed. Only reverse flow in the main-stream direction is 
precluded because of the parabolic marching nature of the solution procedure in this 
direction. 

The pressure can be advanced by rewriting the x-derivative in the continuity 
equation as follows: 



= T p x + T u > 



T 


x 


The p x term can be used to march to the next station. There is no explicit viscous term 
present and any dissipation which may tend to damp pressure oscillations must come 
through a coupling with the other field quantities present in the continuity equation. 
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If equations (5) to (9) are cast in nondimensional form by using free-stream values 
of u . p , and T , and the minimum initial jet width h as the reference length, the 
nondimensional equations become 

for x-momentum, 
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for energy, 
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for continuity, 

|^(pu) + -|(pv) + £<pw> = o 


(13) 


(14) 


for the equation of state, 


p = pT 


and for Sutherland's law, 
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These nondimensional equations are to be solved subject to the following boundary 
conditions: 


y,z — oo: u = 1; v = w = 0; T = 1; p = 1 


= 0: 
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az 
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and initial conditions specifying the particular jet flow problem to be investigated. 

Equations of this type (eqs. (10) to (14)), where the pressure is computed and not 
assumed, have been shown to be singular at M = 1, if the term is treated implicitly, 
and singular at M = 0, if the pressure gradient is calculated in an explicit manner during 
a numerical calculation. (See ref. 6.) If the term is neglected entirely or specified 
from a boundary-layer approximation, then the parabolic march in x can proceed with- 
out difficulty. (See ref. 6.) 

Thus, since the class of problems chosen to test the overall method should avoid any 
of these integration difficulties, if possible, equations (10) to (14) will be solved for a 
supersonic jet issuing into a supersonic free stream. The M = 0 behavior takes place 
near boundaries (where u = 0), and therefore, a free jet problem avoids this singular 
behavior. However, the jet cannot exhaust into an ambient atmosphere since here u = 0 
also. Thus, a jet issuing into a moving free stream is suitable. To avoid any difficulties 
at M = 1, both streams were chosen to be supersonic. 


NUMERICAL PROCEDURE 


Integration Technique 

An implicit numerical procedure was chosen to solve the governing equations for a 
number of reasons. The success of implicit methods on the two- and three-dimensional 
boundary-layer equations implies that they should be efficient for the boundary-layer-like 
parabolic -elliptic Navier-Stokes equations (eqs. (10) to (14)). It is expected that solutions 
may be required at large distances downstream from the initial data plane; consequently, 
large x-steps are desirable. The need to eliminate the step-size restrictions of explicit 
methods leads to a consideration of unconditionally stable methods which are consistent 
in their marching variation. The particular implicit method used in this study is the 
alternating direction implicit (ADI) method of Peaceman and Rachford. (See ref. 5.) The 
ADI method is ideally suited for the solution of equations (10) to (14). There is no stability 
restriction on the step size (see appendix A) and hence, large x-steps are permitted. The 
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method has second-order truncation error in its marching variation, which is also a 
requirement for the type of flow envisioned, since the x-history of the flow must be traced 
accurately at each step. Finally, the method does not require the inversiort of a sparse- 
banded matrix as a fully implicit or Crank-Nicolson method would. Simple tridiagonal 
coefficient matrices which are generated at each step require much less storage and 
time for their inversion in relation to sparse -banded matrices. The method has pre- 
viously been shown to be effective for a set of equations similar to those used in the 
present approach. (See ref. 4.) The use of this technique has also been attempted on 
the parabolic -elliptic Navier-Stokes equations in cylindrical polar coordinates for the 
flow over a cone at incidence. (See ref. 10.) The authors of this work state that an ADI 
procedure was attempted, but failed to work because of difficulties encountered on the 
azimuthal sweep of the ADI method. It is possible that the nonsymmetry of the crossflow 
momentum equations in polar coordinates caused the difficulty. Nevertheless, in Cartesian 
coordinates, equations (10) and (14) are perfectly symmetric in yz and vw, and this 
objection does not arise. 


Coordinate Transformation 


Equations (10) to (14) were cast in finite -difference form and solved for the super- 
sonic free jet problem. Major difficulties were encountered because of the imposition of 
boundary conditions at arbitrary distances into the assumed free stream. The details of 
these results are discussed in reference 13. The exact asymptotic behavior of the solution 
was not known; thus, the conditions in the free stream could not be imposed at any particu- 
lar chosen distance away from the jet orifice since the eventual growth of the jet, as it 
proceeds downstream in the x-direction, will reach this boundary and make the imposed 
conditions incompatible with the flow. To remedy this, a coordinate transformation is 
introduced into equations (10) to (14) which allows conditions at infinity to be imposed 
at some finite distance in the computational domain. This is accomplished by the 
transformation 


V = 


5 = 


Ay 

1 + Ay 
Bz 

1 + Bz 


(15) 


which maps zero onto zero, and infinity in the physical plane onto unity in the plane. 
This linear transformation has been used successfully by other investigators on a variety 
of problems. (See refs. 15 to 17.) More complicated transformations do exist (refs. 18 
and 19), but the transformation (eq. (15)) produced good results; consequently, no others 
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were tried. With equation (15), the governing equations of the problem to be solved 
(eqs. (10) to (14)) are (with £ = x) as follows: 
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These equations are subject to the transformed boundary conditions 
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The ADI procedure is used to difference equations (16) to (20) in x, with the y- and 
z -derivatives expressed as central differences. As an example, by using the notation 

VU = ^Hr + i-”L*-i) (21a) 

(21b) 

the complete differencing scheme for the x-momentum equation is for the first step: 
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where the product of unknown derivative terms which appear in equation (17) have been 
averaged to make the sum of the two ADI steps second order accurate about the point 

H> 

The only difficulty arises from the cross derivatives of velocity present in the y- 
and z -momentum equations. These derivatives cannot be differenced implicitly since the 
tridiagonal structure of the resulting matrices would be destroyed. These derivatives are 
treated in the same manner as all the nonlinear coefficients present in the differenced 
equations. 


Linearization Scheme 

The ADI procedure is a two-step, second-order-accurate, implicit integration tech- 
nique. The first step advances the solution from location (i) to an intermediate stage of 
computation, much like a predictor step; and the second step advances the solution to the 
new location (i + 1), much like a corrector step. The solution generated is then second 
order accurate about the point ^i + (See fig. 1.) This point is not equivalent to the 
intermediate, *, step of the ADI procedure. Hence, a method must be developed to com- 
pute all the nonlinear coefficients (and the cross derivatives) at the i + point. This 
computation can be accomplished by a quadratic extrapolation from the two previous 
x- stations i and i - 1 (ref. 4) by use of the predictor -corrector procedure of Douglas 
and Jones (ref. 20) or by the recursive iteration technique, presented in reference 13, 
common to many numerical boundary-layer procedures. (See refs. 21 and 22.) This 
last method consumes a great deal of time (each iteration is equivalent to a full march- 
ing step), and it was found that the quadratic extrapolation yielded the same results, 
although it required slightly more computer storage. Since computer time was a lim- 
iting factor in some of the calculations to be presented, all the nonlinear coefficients 
and cross derivatives were calculated at the i + w point by extrapolation by using 


Qi 


1+ 2 


j( Ax + ) + ( Ax -)|Qj,k- j( Ax + )Qj,'k 

(Ax.) 




(22) 


where Q refers to any of these quantities, and allowance has been made for variable step 
size in the x-direction. 


Initial and Boundary Conditions 

The solution domain for the free jet problems was the semi-infinite quarter plane 
bounded by the axes of the symmetry of the jet and extending to infinity in both cross 
directions (y,z) in the actual physical plane. The boundary condition difficulties 
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encountered when assuming an arbitrary position for the free stream (ref. 13) were over- 
come by use of the coordinate transformation (eq. (15)). This usage allows imposition of 
the unknown free-stream conditions at the last computational node, that is, at 77 = 1 or 
£ = 1. Conditions on the two lines of symmetry present were fixed by setting the velocity 
normal to the symmetry line equal to zero, and enforcing the symmetry condition of zero 
gradient for all the other variables by fitting a quadratic polynomial through the first three 
grid points, that is, 


4Q i 2 

Q\ . = —h± 

J,i ; 


Q h 


Problems encountered when this relationship was not used have also been discussed in 
reference 15. 

The initial conditions at the data plane representing the jet exit x = 0 were chosen 
in a very rudimentary manner (except in the one case described in the section "Results 
and Discussion” under the heading "Smooth Initial Profile”). Since too many points would 
be necessary to describe the actual merging of the free-stream and duct flows just past a 
jet exit, the initial jet and free-stream u -velocity profiles were represented by two distinct 
flows, separated by a sharp boundary. Computationally, this relationship yields a one 
grid-point discontinuity between uj e ^ and u^. This initial condition is probably the 
most severe that can be imposed while still generating an eventually realistic flow 
description. For simplicity, the initial values of the crossflow velocities in the plane 
x = 0 were set equal to zero. 

The pressure distribution was chosen to be uniform at the free-stream level since 
an unmatched static pressure would undoubtedly produce shocks. These shocks were not 
consciously sought as part of the problem, and the initial conditions were set to try to 
avoid their generation. The streamwise velocity and temperature levels were computed by 
assuming constant total temperature in the jet and free stream and specifying the jet and 
free-stream Mach numbers. 


Solution of Matrix Equation 

By using the extrapolation (eq. (22)) for the nonlinear coefficients and cross deriva- 
tives in equations (16) to (20), the ADI difference approximation for these equations results 
in a linear set of five equations in five unknowns at each step of the ADI marching pro- 
cedure. Rather than make the quite arbitrary choice of the order of solution if a sequential 
technique of solving each equation in turn were elected, the five equations are coupled and 
solved simultaneously. The procedure is quite advantageous because, in addition to 
eliminating the choice of solution order, it models the physics more precisely by allowing 
changes in any variable to be instantaneously sensed by all the others. The resulting 
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matrix equation for the sweep of the ADI procedure where the z -derivatives are differenced 
implicitly has a block tridiagonal structure, one line of which can be represented as 


[ A k]"k,l + [ B k]*k - [ c k]»k-l ' 5 J 
where the unknown vector W k contains the five unknowns 

W k = (u,v,w,T,p) T 


the coefficients 

[A], 

M. 

and [c] 

are matrices, for example, 


a il 


a 12 

a 15 


a 21 


a 22 

a 25 




(23) 


(24a) 


(24b) 


|_ a 51 a 52 a 55j 

and the vector D is the source term in each equation. The components aj^ represent 
the coefficient of the nth unknown from the mth equation at point k and are given in 
appendix B. The inversion of this block matrix is a standard procedure. (See ref. 21.) 


Computer Requirements 

The equations describing the flow were solved in the transformed plane 0 2 77 , £ £ 1. 
A uniform grid spacing in 77 and in £ was used to insure second-order accuracy. 

Two different meshes were used for all the computations, a coarse grid of 21 x 21 points 
and a finer grid of 41 x 41 points. The core storage required for each of these programs 
was 60Kg and 15 1 X 3 , respectively. The time required at each x station to process 
one unknown at one ( 17 ,?) point was approximately 0.00165 second per unknown per point 
on a CDC 6600 computer. This time is. slower than the value of 0.0005 for a sequential 
ADI scheme used on a time asymptotic Navier -Stokes solution (ref. 23), but compares very 
well with the approximate value of 0.0313 computed from the 30 seconds of CDC 7600 com- 
puter time, quoted in reference 10, for a 50 X 23 grid (assuming a CDC 7600 computer is 
6 times faster than a CDC 6600 computer). Thus, to compute one x-step of the PENS 
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procedure took 3.6 seconds for the 21x21 grid and 13.9 seconds for the 41 x 41 grid. 
Typical runs consisted of approximately 100 steps. 

RESULTS AND DISCUSSION 

The supersonic jet flow problem was solved for a number of different free-stream, 
geometrical, and computational conditions; these are summarized in table I. In all cases, 
(1) the Prandtl number Np r was set equal to unity, (2) the reference length was set equal 
to the minimum initial orifice width of the lxl square jet, (3) the pressure distribution 
was taken to be uniform at the free-stream value throughout the domain at x = 0, (4) the 
initial velocity distribution was discontinuous except in the case labeled hyperbolic tangent 
profile in table I, (5) the total temperature in the entire field was set equal to a constant 
at x = 0, and (6) the ratio of initial jet to free-stream Mach numbers was always 1.5. 

The marching procedure followed for all runs was to initially use a step size Ax 
which was small compared with the core length of the jet being calculated, and then incre- 
ment the x-step size occasionally in the march in order to be able to reach significant 
downstream distances from the jet orifice. 

The main purpose of the test cases was to demonstrate the applicability of the 
equations, numerical procedure, and computer code. For axisymmetric or two- 
dimensional flow, the usual boundary-layer assumptions hold, and a standard free shear 
layer computation can be used to generate the flow. These two limiting geometric cases 
(two-dimensional and axisymmetric) are used as references for all the three-dimensional 
calculations. A square jet and increasing aspect ratio rectangular jets were computed. 

The square jet will be seen to behave much like an axisymmetric jet and the rectangular 
jets reproduce two-dimensional behavior for a part of the distance away from the orifice, 
but ultimately all the jets resemble axisymmetric ones at large distances downstream. A 
circular jet was also approximated in the three-dimensional Cartesian grid, and the dif- 
ferences between this and the previously computed flows are noted in addition to a com- 
parison with the free shear layer calculation. A final case of interacting jets is shown to 
demonstrate the range of problems which can be handled. 

Thus, in each case (except the interacting jets) a comparison will be made between 
the PENS results and the results obtained by using a classical shear layer jet (SLJ) code 
(ref. 24) of high accuracy. 

In all the cases to be presented since Np r = 1 and the initial total temperature 
was constant, it was expected that this constancy would remain throughout the calculation. 
During the initial stages of jet development, there was about a 5 -percent deviation from 
this constant, but when the calculation began to approximate the true asymptotic boundary- 
layer behavior, described in the subsequent sections, the constant total temperature relation 
was well approximated. 
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Square and Axisymmetric Jets 

The variation of the velocity on the center line of the computed Mj = 1.8, Nj^ e = 10 
jet with distance downstream of the initial orifice location is shown in figure 2. This 
gross feature of the flow is an excellent indicator of both the qualitative and quantitative 
aspects of the calculation method. On the low storage 21x21 grid, two transformations 
corresponding to different values of the constant A and B in equation (15) were used. 
The relation between points in the physical and computational plane is shown in table II. 
Thus for A = B = 4, only 5 points are originally in the jet at x = 0 compared with 
11 points in the jet when A = B = 2. The center-line velocity is noticeably affected as 
seen in figure 2, with the A = B = 2 case being shifted toward the shear layer result. An 
important detail which cannot be seen on the scale plotted in figure 2(a) is an overshoot of 
the center-line velocity in all cases when A = B = - was used with this grid. Although 
the amount of overshoot was not particularly high (0.3 percent), it did indicate the presence 
of an inaccuracy. When the constant in the transformation was changed to A = B = 2, 
a monotonic decay of the center-line velocity always occurred. 

For even greater accuracy, a 41 x 41 grid was used for the square jet, along with a 
transformation of A = B = 2. Again, a noticeable shift in the center-line decay curve 
toward the axisymmetric value is apparent. In an attempt to reproduce the SLJ code 
results used as a comparison in all these cases, a pseudo-circular geometry was con- 
structed on the 41 x 41 grid. The resulting "circle" is shown in figure 3. The difference 
in area of the "circle" from a true circle is about 1 percent. The resulting center-line 
decay was indeed close to the SLJ value. (See fig. 2.) 

One reason for the discrepancy of the approximate circle from the SLJ case and, 
in fact, a reason for the continual approach of the four curves ^2 1x21, A = B = ^; 

21 x 21, A = B = 2; 41 x 41, A = B = 2; "circle") to the SLJ case is the difference in 
initial momenta (or for these cases, area) of each jet at x = 0. Since a discontinuous 
profile was assumed for the initial velocities, there is an additional amount of momentum 
contained between the last grid point at the jet edge, and the next grid point which is in the 
free stream. This additional momentum delays the velocity decay somewhat, as seen by 
the displacement of the curves in figures 2(a) and 2(b). As the grid size is decreased 
either by adjusting the transformation constants or by increasing the number of nodes, the 
initial momentum contained within the one grid step beyond the edge of the jet decreases 
and gives the results shown. The significant difference between the square and approxi- 
mate circle curves is caused by the 22-percent reduction in area from a square to its 
inscribed circle. However, there is again an incremental area difference between the 
approximate circle and the axisymmetric case used for comparison. Since the SLJ code 
had 256 grid points, its approximation to a circle was very good. On the other hand, the ' 
PENS code had about 8 percent more momentum beyond the jet than was contained in it. 
This could certainly account for the shift of the PENS curve away from the SLJ curve. 
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As a final check on the ability of the code to compute correctly, at least for global 
features of the flow such as center-line decay, a comparison was made at Mj = 1.8 of 
square jets with Reynolds numbers of 10 and 1000. This comparison is shown in figure 4. 
The curves are identical to the scale used for the plotting; table in presents the computed 
values, and the difference is seen to be usually in the third significant figure. Hence the 
computation reproduces the correct Reynolds number scaling for an axisymmetric jet 
(See ref. 2.) 

A more severe test of the accuracy of the computer code is a comparison of actual 
velocity profiles rather than an integrated feature of the flow, such as the center-line 
decay rate. Figure 5 shows a comparison of the u velocity profile calculated for a cir- 
cular jet with a Reynolds number of 10 and a jet Mach number of 1.8. Since it has already 
been seen that the PENS deck initially carries more momentum and will not exactly repro- 
duce the x-wise behavior of the SLJ deck, two x-positions corresponding to approximately 
the same center-line velocity were compared. The results of the SLJ curve are for 
£ = 2.3, whereas the ^ = 2.5 points are plotted for the PENS curve. Since these 
x-positions are in the asymptotic similarity region of the jet, equal center-line velocities 
should produce comparable profiles; figure 5 verifies this assertion. The two profiles 
almost coincide. It should also be noted that all the y-points used in the PENS calculation 
which fit within the scale of the plot are shown, but for the highly accurate SLJ calculation 
175 points were within the plotting scale; thus only a few are shown. 

The most sensitive of all the computed variables is the velocity component normal 
to the free stream in the cross planes. In the PENS calculation where the pressure gradi- 
ents in the cross planes are actually computed and act as a driving force for v and w, 
these crossflows are highly influenced by any irregularities in the flow field. This is an 
important feature of the parabolic-elliptic Navier-Stokes calculation procedure and points 
out a significant difference between this method and standard shear layer calculations. 

The boundary-layer assumptions give the imposed pressure through the shear layer from 
the reduced normal momentum equation. The normal velocity is then computed from the 
continuity equation used in an auxiliary role to conserve mass. In the PENS procedure the 
momentum equations are used to compute the velocities and continuity gives the pressure. 
The roles of the equations and the variables they generate are completely reversed in the 
two techniques; see reference 25 for a further discussion. 

The computed crossflows for the Nr 6 = 10, Mj = 1.8 circular jet are shown in 
figure 6. Here the agreement between PENS and SLJ is still good in the major portion of 
the jet, but deteriorates in the outer region. Also, there are obvious ripples in the PENS 
solution. This behavior is caused by oscillations in the computed pressure around the 
initial orifice location, ^ = -i, and can be traced back to the initial discontinuity in the 
velocity and temperature profiles. An initial fluctuation in the pressure develops 


20 



(ref. 13) but eventually decays so that pressure in the cross plane is effectively a constant 
(the free-stream value). However, superimposed on this constant value are oscillations 
which are essentially undamped. (See appendix A.) Since the cross velocities are small 
(« 1CT 3 ), relatively small oscillations in the pressure will be felt through the pressure 
gradient terms in the cross momentum equations. However, the results are quite good 
despite the oscillations present. (See fig. 6.) 

The effect of the pressure terms can be seen more clearly if a comparison of veloc- 
ity profiles is made further upstream in the jet flow near the end of the core region. 

Figure 7 shows the u-velocity for both PENS and SLJ calculations at x = 0.13 when the 
center -line velocity ratio (u^ - U^y^u^ - U M ) is still equal to one. Again, the two 
curves are almost coincident. However, the comparison of v is striking. (See fig. 8.) 

The SLJ calculation shows the typical jet crossflow profile. The PENS result is very 
different. Here the v-profile almost precisely follows the cross-plane pressure gradient 
set up by the discontinuous initial profiles. Obviously, the pressure and cross -velocity 
effects combine to give the correct behavior of the main flow as evidenced by the agree- 
ment in figures 5 and 7 and the eventual agreement of the cross velocity further down- 
stream (fig. 6). This demonstrates that the PENS and SLJ procedures react differently to 
the initial profiles. In the SLJ code all discrepancies of the flow are absorbed by the 
normal velocity through the continuity equation. The PENS code has an additional degree 
of freedom (the calculation of the pressure) and so initially adjusts differently, but even- 
tually computes a standard shear layer flow. 

Additional calculations were made for Mj = 1.8, Nr 6 = 1000, and Mj = 7.5, 

Nj^ e = 1000. The center-line velocity decay for the Mj = 1.8, Nr 6 = 1000 jet is shown 
in figure 9. The 41 X 41 square jet has already been compared for a Reynolds number 
scaling with the Nr 6 = 10 jet in figure 4. The center-line velocity plot for the Mj = 7.5, 
Nr 6 = 1000 case is shown in figure 10, the results being similar to the others thus far pre- 
sented. A comparison of the PENS and SLJ u-velocity profiles is given in figure 11 for the 
Mj = 1.8, Nj^ e = 1000 circular case. Again, these comparisons show good agreement 
between PENS and SLJ calculations. However, no comparisons can be made of the cross 
velocity due to the behavior of the calculated pressure. Whereas the Nr 6 = 10 jet 
exhibited oscillations only in the YZ plane, for the Nr 6 = 1000 cases, an oscillation 
also appears in X. It is small (usually third or fourth decimal place) but significant. On 
alternate x-steps the pressure oscillates about the free-stream value of 1, and causes the 
y- and z-pressure gradients to change; hence, v and w change sign at successive 
x-steps and give nonphysical results. It is concluded that at this higher Reynolds num- 
ber, the initial discontinuity affects the solution severely, much more than it did in the 
N Re = 10 case. 
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Smooth Initial Profile 


To test the solution procedure for less drastic initial profiles than the discontin- 
uous ones used previously, a smooth, hyperbolic tangent initial profile was used for u in 
the PENS code. Again, the total temperature was assumed to be constant at x = 0, and 
the other variables were initialized as before. The initial u profile had no core region 
near the axis and so diffused rapidly once the PENS marching procedure was initiated. 

The center-line decay rate is shown in figure 12 for the Njj e = 10, Mj = 1.8 PENS 
approximate -circular and axisymmetric SLJ cases. The agreement is very good. The 
u-velocity profiles compare very well with each other over the entire x-range computed. 
The normal velocity components behave in the same way as before, except that there are 
no oscillations in the solution, unlike those shown previously. 

Even more important than the elimination of the low Reynolds number oscillations, 
the smooth initial profile allows the calculation of sensible crossflows at the higher 
Reynolds number of 1000. (See fig. 13 for a comparison of PENS and SLJ calculations.) 
Using the previous discontinuous profile produced meaningless values of v and w 
because of ripples in the calculated pressure. These ripples, although not completely 
eliminated, are reduced by more than three orders of magnitude when the smooth initial 
data were used, and their effect on the crossflow is greatly diminished. 

Varying Aspect Ratio Rectangular Jets 

The three-dimensional capability of the PENS code can be tested by varying the 
aspect ratio of the initial profiles, from the square one used previously, to rectangular 
profiles of large width-height ratios. These high-aspect-ratio jets should reproduce two- 
dimensional behavior near the center line for some distance downstream of the initial 
station until the effect of the initial edge of the jet propagates in far enough to change the 
flow to a three-dimensional one. 

An excellent way to display this phenomenon is, again, to plot the streamwise veloc- 
ity on the center line. Figure 14 shows, for Nj^ e = 10, Mj = 1.8 jets, the center-line 
velocity decay. A 4 x 1 rectangular jet and an 8 x 1 jet are displayed, along with a shear 
layer calculation of true two-dimensional jet. The two rectangular jets closely follow the 
SLJ results far downstream before three-dimensional effects alter the decay curve and 
finally cause it to decay in an axisymmetric manner. Initially, the two rectangular jets 
follow exactly the same curve. Beyond some location ^approximately ^ = 3 in fig. 14^, 
the outer edge of the initial profile is sensed earlier in the 4x1 case than in the 8x1 
case, and the center line accordingly decays more rapidly. In the final stage of their 
decay, both jets behave axisymmetrically. In fact, they very nearly simulate the expected 
behavior with regard to their initial momenta. Since the static temperature at this stage 
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is very nearly equal to the free-stream value throughout the jet, the incompressible decay 
rate for axisymmetric jets, 


Uoo 87 r x/h 

(see ref. 2, eq. (10.19)) where K is the nondimensional initial jet momentum, can be 
used as a comparison. This relation shows that for equal ratios of initial momenta to 
x/h, the center-line velocity should be the same. Since the momentum in the 8 x 1 jet is 
twice that of the 4 x 1 jet, the velocity of the 8 x 1 jet should be equal to that of the 4 x 1 jet 
at half the x/h of the 8x1 jet. The values shown in figure 14 very nearly confirm 
this ratio. 

One further point about the results of figure 14 should be noted. The computed two- 
dimensional decay rate is not the usual value of x - -*-^ of two-dimensional incompressible 
flows. (See ref. 2.) This was noted previously by Pai (ref. 26) who attributed it to the 
presence of heat transfer. 

The u- and v-velocity profiles computed by the PENS method for an 8 x 1 rectangle 
(at x/h = 1.5) and the SLJ method (at x/h = 1.3) are shown in figures 15 and 16. The 
different x-locations were chosen to equate the center-line velocities as previously 
described. Close to the axis, both calculations agree very well, but in the outer region of 
the jet, there is a substantial difference. The PENS jet has spread further into the free 
stream. 

As these x-locations are not in the asymptotic decay region, exact comparisons 
should not be expected and the excess spread of the PENS jet, calculated from a true 
three-dimensional procedure, can be explained. To reach the given value of center-line 
velocity shown, the PENS calculation has to proceed further downstream and hence will 
spread somewhat more than the SLJ jet. 

The true three-dimensional nature of the PENS calculation also contributes to the 
additional spread by removing the constraint of only one direction from which entrained 
momentum may be drawn. The effect of the second cross stream direction is always felt, 
although initially it is small, and allows another degree of freedom by which the jet may 
spread. 

For the higher Reynolds number jets, the same quantitative behavior was obtained 
in the rectangular cases as in the square and circular cases. 

Interacting Jets 

As an example of the wider range of geometries which can be handled by the PENS 
calculation method, a case of the flow generated by two interacting jets was computed. 
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Two square jets of width h, separated by a distance h, were computed under the same 
flow conditions as previously used in some of the PENS test cases, that is, M*, = 5.0, 

Mj = 7.5, N Re = 1000. (See table I.) 

This sort of geometrical configuration is of interest in the aerodynamic integration 
of exhaust nozzles on the current F-15 configuration into the fuselage. This integration 
entails the present configuration containing two axisymmetric nozzles altered as shown 
in figure 17 to contain two two-dimensional (rectangular) nozzles. No attempt was made 
in the calculation to duplicate the actual conditions of an F-15 nozzle; one of the cases 
calculated previously was simply run with altered initial geometry to simulate the inter- 
acting jets. No difficulties were encountered with this computation. 

Figure 18 shows the behavior of the streamwise velocity on the symmetry line 
between the two jet exits, the square and 4x1 rectangular jet results being included for 
reference. Initially, on the center line the velocity is the free-stream value and initially 
increases as the influence of the jet is felt. Eventually though, the offset jets merge and 
the center-line velocity starts to level off and then decay quickly, axisymmetric behavior 
being assumed. Figure 19 shows the computed streamwise velocity profiles on the sym- 
metry line, and figure 20 shows these computed velocities in a perspective plot for the 
half plane y ^ 0 of the flow domain. These results are qualitatively correct; however, 
no. comparisons can be made with shear layer calculations since this three-dimensional 
flow is beyond the scope of boundary-layer theory. 

Possibility of Extending the PENS Method to Turbulent Flows 

All the comparisons of PENS and SLJ calculations which have been made thus far 
have been for laminar flows. It will now be necessary to include the appropriate turbu- 
lence model for these flows in order to extend the applicability of the procedure. A 
second- order turbulence closure method has already been used with the incompressible 
PENS approximation in the calculation of submarine wakes without reported difficulty. 

(See ref. 12.) Hence, in principle, a compressible model should pose no insurmountable 
problems, although the task is by no means simple. In fact, the inclusion of turbulence 
modeling into the solution procedure may very well take better advantage of the present 
PENS procedure since the coupled solution technique is almost certainly necessary in the 
turbulent case (ref. 27), but is probably not essential for the solution of the laminar flow. 
Although the models complicate the governing equations, perhaps making them more 
difficult to solve, the effective Reynolds number of the flow is reduced by the addition of 
turbulence; and the lower Reynolds number cases which have been computed yielded 
significantly smoother results. 
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CONCLUDING REMARKS 


In all cases where comparisons could be made with shear layer calculations, the 
parabolic-elliptic Navier-Stokes calculations produced results which compared well, not 
only with gross characteristics of the flow, such as center-line decay rate, but also with 
detailed profile calculations, except in the initial stages of the flow development where the 
procedure calculates a pressure gradient across the shear layer, to which the crossflow 
velocity reacts differently than a standard shear layer flow. 

The alternating direction implicit method has been successfully applied in a numeri- 
cal solution of these PENS equations. A coordinate transformation was introduced which 
allowed the imposition of free-stream (that is, infinity) boundary conditions at a finite 
location in the computational domain. This transformation eliminated any difficulties 
caused by an arbitrary choice of free-stream location in the untransformed plane but 
created some difficulties of its own. Too large a coordinate stretching in a region of 
significant flow proved to yield inaccurate and nonphysical results. 

Discontinuous initial conditions were used and produced smooth results at low 
Reynolds numbers, but in the higher Reynolds number cases oscillations in the pressure 
caused a deterioration of the calculated crossflow. This problem was alleviated by using 
a smooth initial profile. Thus, if experimental data were used at the initial station, it is 
assumed that no difficulties would ensue, since the data would presumably be smooth. 

Computed results have been presented for various jet configurations. Although the 
character of the computed flow can be predominantly axisymmetric in the case of the 
square jet or approach a two-dimensional limit as is the case for the 8x1 rectangular 
jet, no prior assumptions to this effect are required. The solution computes the particu- 
lar flow under investigation as though it were three-dimensional and allows the initial 
geometry of the prescribed jet to determine the ultimate nature of the solution. 

A final configuration studied consisted of two square jets symmetrically disposed 
about a symmetry axis. This is a simple geometric model of the exhaust of a proposed 
F-15 modification and demonstrates the range of geometries which can be handled. 

The parabolic -elliptic Navier-Stokes equations have been shown to model adequately 
flows where the two crossflow directions are of equal importance in a predominantly con- 
vective viscous flow. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
August 12, 1976 
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APPENDIX A 


STABILITY AND DAMPING OF THE ADI METHOD 


Any of the governing equations used to solve the jet problem (eqs. (16) to (20)) can 
be linearized to give the simple form 


9Q 

dx 


+ 


a^ + b^ = 
dy dz 




(Al) 


Note that the following analysis is only for one equation in one unknown, and not for the 
entire coupled system (eqs. (16) to (20)) for which the analysis is much more difficult. 
However, as each equation can be thought of as influencing one particular variable more 
than any other, that is, momentum equations and corresponding velocities, energy equation 
and T, and continuity equation and p, some insight into the numerical behavior of the 
ADI method can be obtained. 

By using the ADI method, the difference form of equation (Al) becomes 


%,m ^i,m 


“ , + a6 v Q, + b6 7 Q, 

(Ax/2) r%m z^,m 


Qz,m Qz,m . „ R ^n+1 , 

’-r- 7 — - 2 — + ao v Q, „ + bo 7 Q 7 __ 

(Ax/2) y t,m z l,m 




(A2) 

(A3) 


where l and m are the discretized y- and z-positions and *, n, and n+1 are the 
discretized x-positions. Performing a Von Neumann stability analysis using 


Q 


l, m 


= Q n e*(' 


Zky Ay+mk z Azj 


in equations (A2) and (A3) gives 


Q* - Q n + ic y Q n sin 0 y + ic z Q* sin <j> z = jS y Q n (cos 0 y - 1) + /3 z Q*(cos cp z - 1) 

Q n+1 - Q* + ic y Q n+1 sin + ic z Q* sin $ z = /3yQ n+1 (cos 0 y - 1) + 0 z Q*(cos <p z - 1) 
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where 


and 


, _ aAx 
y 2 Ay 




a Ax 
Ay 2 


0 y = k y*y 


_ aAx 
z 2Az 




aAx 

Az 2 


<t> z = k z Az 


Therefore, 

G * _ Q* _ 1 ‘ ic y sin ^y " %( 1 “ cos ^y) 

Q n 1 + ic z sin 0 Z + / 3 Z (1 - cos <f > z ) 


G n+1 = Q n+1 = 1 - ic z sin 0 g - fl z (l - cos <ft z ) 

Q* 1 + ic y sin (p y + /3 y (l - cos <£ y ) 

and therefore the amplification factor G for the complete step is 

G = G*G n+1 = 1 - ic y sin <py - M 1 ~ cos fy) 1 ~ ic z sin 4> z - M 1 - cos flz) 
1 + iCy sin <py + /3y^l - cos 1 + ic z sin <p z + /3 Z (1 - cos $ z ^ 


(A4) 

Hence |G| = 1 and the method is stable for the convection-diffusion equation (Al). 

If equation (Al) is simplified so that there are only convective terms, that is, a = 0, 
then it becomes a simple linearized model of the continuity equation (where Q = p). For 
this case, the amplification factor (eq. (A4)) becomes 


G = 


1 - iCy sin cpy 1 - ic z sin $ z 
1 + iCy sin (py 1 + ic z sin 4> z 
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■The magnitude of G is therefore 

l G | = 1 

exactly, and the numerical scheme introduces no damping. Oscillations which may occur 
in the solution do not decay and persist undamped (unless affected by a coupling to the 
other equations where viscous terms are present). This is not the case for any other 
implicit scheme except a Crank-Nicolson differencing, which is equivalent to the ADI 
method within the second-order truncation error. (See ref. 28.) 

A general implicit scheme (ref. 28) for equation (Al) with a = 0 gives for the 
amplification factor 

1 - i(l - 9) (Cy sin </> y + c z sin <£ z ) 

1 + i0(c y sin <py + c z sin 0 z j 

gives Crank-Nicolson differencing and 9=1 gives a fully implicit dif- 
The magnitude of G is 

U[l- 8) 2 (c y sin 4>y + C z Sin 4> z f (A5) 

1 + # 2 (Cy sin 0y + c z sin <£ z ) 2 

To find the damping, look at the long wavelength disturbances by letting <py and <p z 
approach zero simultaneously, and expand the denominator of equation (A5) to get 

| G | = 1 + (1 - 9) (c-y'fiy + ^z^z) " 1 ” ®^( c y^V ^z^z) + • • • 

Thus 

|g|= 1 + (1 - 20)( Cy 0y + C Z 0 Z ) 2 + . . . 

and for only 9 = A, that is, Crank-Nicolson differencing, is the amplification factor iden- 
tically equal to unity. For 9 = 1, a fully implicit scheme, significant second-order damp- 
ing occurs. 


where 9 = 
ferencing. 
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'i 

! 

| 

MATRIX COEFFICIENTS OF THE LINEARIZED DIFFERENCE EQUATIONS 

In equation (23), for the first half of the ADI procedure, from a location i to the 
intermediate * location, the components of the matrices {Aj, [b], (CJ, and the vector D 
are 

a ll = + 2B ^ " S)*y%e B(1 - ?) 2 / a ? - B 2 (l - O 4 ^/ N Re 

a 12 = 0 
a 13 =0 

a 14 = -| B2 d-?) 4 |^||/ N Re]/ A ? 
a 15 = 0 
a 21 = 0 
a 22 = a ll 

a 23 = | AB(1 - „> 2 (1 - ?) 2 Hlf/Npely^ 

a 24 = -{{[b 2 * 1 - « 4 |f + AB d - ?) 2 d - 1> 2 |fjw/ N Ite| A 
a 25 = 0 
a 31 = 0 

a 32 - -|j AB <1 - 0) 2 d - « 2 lrl^/ N Ee^5 

a 33 - Ijp* + ! B(1 - C)</nJ b(1 - f) 2 Ac - | B 2 d - r} 4 ip./^: 2 * M |t/a k 
a 34 ■ AB <1 - n)\ 1 - a 2 §* - B 2 d - C) 4 4 

a 35= [iBd-oywJA 
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>41 “ - [e~MB 2 <1 - ?)* |f/Nite] A? 

>42 = - ^ - » 4 |2 + AB(1 - „) 2 <1 - C ) 2 A/h 

>43 - E^ixpd - C ) 4 9w - i AB(1 - ») 2 (1 - C ) 2 f|^A 

>44 ■ |[pw + 2B(1 - ?)B/N Re ^ Pr ]B(i - 0 2 /a? - B 2 u - 0 4 [p/a<; 2 + 1 It 
a 45=-|^ B(1 - C) V^ 


>51= 0 
a 52 = 0 

>53 = I B <> - « V* 


a 54 = - i B(1 - C) 2 pw/a< 

>56 * j B(1 - !) ! »A 

b n = 2pu/A{ + + |2B 2 (1 - S^M/NRejA 2 



•>15 = ( 2 /*4)/>' M « 

» 21 *° 
b 22 = b ll 
b 23 = 0 
b 24 = 0 
b 25 = ° 
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b 31 = 
b 32 = 
b 33 = 
b 34 = 
b 35 = 
b 41 = 
b 42 = 

b 43 = 

b,, = 
44 

b 45 = 
b 51 = 
b 52 = 
b 53 = 
b 54 = 
b 55 = 
c ll = 
c 12 = 
c 13 = 
c 14 = 
c 15 = 


0 

0 



0 

0 



0 


-2pu/A| + 

2u/a£ + 

m/a ? 2 


•||pw + 2B(1 - OM/N R e]B(l - ?) 2 /a? - B(1 - ?)' 


0 


0 


" a 14 


0 


1 djj. dT 
4 9T 8? 
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c 21 = 0 


c 22 " C 11 


c 23 “ " a 23 


c 24 " _a 24 


c 25 “ 0 


C 31 _ 0 


c 32 _ " a 32 


c 33 _ 


c 34 _ 


f[pw + §B(l- ?)p/n r J]b(1 - 0 2 /aC - | B 2 (l - C) 4 4 /i/a^ 2 - 


-a 
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c 35 “ ” a 35 


C 41 " “ a 41 


c 42 _ -a 42 
c 43 = " a 43 


c 44 = " - 


pw + 2B(1 - Om/n- 


Re/ N Pr| 


B(1 - 0 2 Ac - te*(i - ?)‘ 1 |m/a^ 


la 

9T 


9T 

w 


: 45 = 

" a 45 

51 = 

0 

52 = 

0 

53 = 

" a 53 

54 = 

" a 54 

55 = 

" a 55 


32 



APPENDIX B 


d i * 2pu u S *h+ + ^2 


Pj )k y^4+ - pv + 2A(1 - 77 )p/N Re A (1 - 7 7 ) 


3 u 
97] . 




9r? 2 2 9T197, ar,| j k 87 , 87 ,^] 


2pu v j,k/* + ' A(1 ' ^ ^2 2 j . - [ PV + f A(1 ' ^Re] A(1 " ^ f [ 

' iVA oo j - J J J 


■ A 2 n «\ 4 4 .. a 2 v , 2 9u far avl 1 , 9v 9T 1 ,, /N 

+ a (l - 77 ) — m —~n + o-^ferid + -sr^r 11 /N 


3 d7]2 j,k 3 8T H a 7 ? |j,k dTld71 \j,: 




3 ^9r, 9? T 9TI2 9? 9 t) j)k 3 3? drj].^ 


d 3 = 2pu - pv + 2A(1 - T))p/N Re A(1 - t ?) 2 ^ 


+ A 2 (l- 77 )%^ + f*|T 1 A; 


^2 2 9T197, 977 j>k 


AO /1 Ai 1 9 2 v . 9p/l 9v 9T 1 1 9T 9v 1 \ 

AB ( 1 - t?) (1 - ?) 3 ^ d £ + aT ( 2 9? 977 . . " 3 9? 977 

\ ^ > K j > k / 


d 4 = 2pu Tj >k A| + - [pv + 2A(1 - C)p/N Re /N pr J A(1 - 77 ) f± 


^2 3TI 1 


(y - 1 ) . ,, .2 9p „ i / A , .2,, A 9 2 T 9p 3T 3T 1 

+ — A(1 - 77 ) v -t=- - 2up^ / A S + + A (1 - 77 ) / p — s- ™t-v 

y 9t ? 4 t ],k/ + ] 37,2 9T 9?7 977 

J > K J i .k J > K 


N Re / N Pr 


, P *2/., \4 /9u 9u 1 4 9v 9v 1 9w 9w 1 \ 

+ E ^ A(1 -’ 1) feij *JKi 5 ; jjk + ir 5 ?J 


2 , ^ 2 /av awl* 2 3w 9v|^ 


+ ABU - u)“(l - / N Re 

\9? 8 1 ; , 3 9t ? i lr / Ke 
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d 5 = 


= 2 (> 


UP J,k 


pu T i,k * 


pu: 


h* 


'H + - A(i - n)‘ 


dtj 


j 5 k 


-pv|X | 
817 1 


i,k 


8v 

+ P 8? 


1 

j,k| 


where the unsubscripted field variables are evaluated at the point (i + j, k\ by the extra- 
polation (eq. (22)). The second half of the ADI step, from location' * to i + 1, has ele- 
ments similar to these, but with the roles of 77 and £ reversed. 
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TABLE L - VARIOUS JET CONFIGURATIONS COMPUTED 


! N Re = 1°’ 

N Re = 1000; 


Configuration 

1 M = 1.2; 

M 00 = 1.2; 




Mi = 1.8 



Square jet; 21 x 21 grid; 

A-B-J 

Square jet; 21 x 21 grid; 

A = B = 2 

Square jet; 41 x 41 grid; 

A = B = 2 

4x1 rectangular jet; 41 x 41 grid; 

A = B = 2 

8x1 rectangular jet; 41 x 41 grid; 

A = 0.52; B = 2 
Circular jet; 41 x 41 grid; 

A = B = 2 

Circular jet; 41 x 41 grid; 

hyperbolic tangent profile 
Interacting square jets; 41 x 41 grid; 
A = B = 2 




TABLE H.- EFFECT OF TRANSFORMATION CONSTANT ON LOCATION 


OF GRID POINTS IN THE PHYSICAL PLANE 
FOR A 21 x 21 GRID 




crlX 


TABLE m.- CENTER-LINE VELOCITY RATIOS FOR Mj = 1.8 SQUARE JETS 
FOR TWO DIFFERENT REYNOLDS NUMBERS 



0.000 

.015 

.02 

.025 

.03 

.04 

.05 

.06 

.08 

.10 

.15 

.20 

.30 

.40 

.60 


“t( N Re = 10 ) 

"t( N Re ‘ 100 °) 

1.0 

1.0 

.9956 

.9952 

.9846 

.9840 

.9665 

.9664 

.9432 

.9437 

.8887 

.8866 

.8318 

.8276 

.7764 

.7689 

.6761 

.6644 

.5918 

.5878 

.4418 

.4398 

.3485 

.3499 

.2427 

.2447 

.1852 

.1874 

.1252 

.1273 
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.1 .3 1.0 3.0 10. 

x/h 

Figure 2.- Streamwise center-line velocity decay for square and axisymmetric jets. 
N Re = 10; Moo = 1.2; Mj = 1.8. Flags denote superimposed values. 
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1.000 1.025 1.050 1.075 

u/u 

7 oo 

Figure 5.- Streamwise velocity profiles for PENS and SLJ calculations. 


Nr 6 = 10; M m = 1.2; and Mj = 1.8. 
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v / u « 

Figure 6.- Normal velocity profiles for PENS and SLJ calculations. 
N Rp = 10; M*, = 1.2; and Mi = 1.8. 








o 21 x 21 grid; A = 2; square jet 
□ 41 x 41 grid; A = 2; square jet 



x/h 

Figure 9.- Streamwise center-line velocity decay for square and axisymmetric jets. N^ e = 1000; 
Moo = t.2; and Mj = 1.8. Flagged symbol denotes superimposed values. 




Figure 10.- Streamwise center-line velocity decay for square and axisymmetric jets. 

N Dq = 1000; Moo = 5.0; and M, = 7.5. 



Figure 11.- Streamwise velocity profiles for PENS and SLJ calculations. 
N Re = 1000 ' M °° = 1.2; and Mj = 1.8. 






Figure 12.- Streamwise center-line velocity decay for axisymmetric jets started with smooth 
initial profiles. N^ e = 10; M m = 1.2; and Mj = 1.8. 
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Figure 13.- Typical normal velocity profile in axisymmetric jet 
calculated using PENS procedure with smooth initial data. 
N r = 1000; = 1.2; and = 1.8. 
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□ Two-dimensional shear-layer calculation; ^ = 1 
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Figure 15.- Streamwise velocity profiles for PENS and SLJ calculations. 
N Re = 10 * M °° = *.2; 311(1 M j = 1,8> 





Figure 17.- Current F-15 configuration and proposed modification with two-dimensional nozzles. 
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Figure 18.- Streamwise center-line velocity decay for interacting jets. 
Square and rectangular jets shown for comparison. 
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Figure 20.- Perspective plots of interacting jet streamwise velocity at various x-stations. 
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